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ABSTRACT 

The median observed velocity width vgo of low-ionization species in damped Lya systems is close 
to 90kms~^, with ~ 10% of all systems showing vgo > 210kms~^ at 2: = 3. We show that a relative 
shortage of such high-velocity neutral gas absorbers in state-of-the-art galaxy formation models is a 
fundamental problem, present both in grid-based and particle-based numerical simulations. Using a 
series of numerical simulations of varying resolution and box size to cover a wide range of halo masses, 
we demonstrate that energy from gravitational infall alone is insufficient to produce the velocity 
dispersion observed in damped Lya systems, nor does this dispersion arise from an implementation of 
star formation and feedback in our highest resolution (~ 45 pc) models, if we do not put any galactic 
winds into our models by hand. We argue that these numerical experiments highlight the need to 
separate dynamics of different components of the multiphase interstellar medium at z = 3. 
Subject headings: galaxies: formation — galaxies: kinematics and dynamics — intergalactic medium 
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1. INTRODUCTION 

Damped Ly-alpha absorbers (DLAs) provide us with 
a high-resolution probe of galaxy formation processes up 
to redshift z ~ 5, allowing us to study the neutral hydro- 
gen distribution and kinematics of young galaxies and 
their surroundings on scales from several pc to several 
kpc, although detailed information about the spatial ex- 
tent of the absorbing gas cannot be readily extracted 
from the observed velocity line profiles. DLAs contain a 
large fraction of high-redshift baryons potentially avail- 
able for star formation (SF; Prochaska et al. 2005), how- 
ever, their in-situ SF efficiency appears to be a factor 
of 20-30 lower than in the present-day galaxies of com- 
parable gas surface density (Wolfe & Chen 2006). Even 
though DLAs are presumably associated with fairly com- 
pact structures, they cover about a third of the sky, so 
there is a high probability that a line of sight to a re- 
mote quasar will contain such an absorber. Currently, 
just over a 1000 such systems have been studied, and 
their number will undoubtedly increase in coming years 
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One of the unsolved long-standing puzzles in our un- 
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derstanding of DLAs is their large neutral gas velocity 
dispersion (Prochaska & Wolfe 1997). The median value 
of Vgo, the velocity interval encompassing 90% of the op- 
tical depth, in absorption lines of low ions associated with 
cold neutral gas (such as Sill, NI, Fell, etc.) is close to 
90kms~^, and 10% of all systems have vgo > 210 km/5. 
Observations indicate that there is virtually no correla- 
tion between the column density Nm of the absorber and 
its absorption line velocity width, with even lower column 
density systems at the DLA threshold A^hi = 10^^*^ cm~^ 
demonstrating velocity widths up to 400kms~^. On the 
other hand, a correlation is observed between the veloc- 
ity widths and metallicities of the DLAs. In fact, the 
correlation is remarkably tight between the equivalent 
width of the Si II 1526 transition (a kinematic diag- 
nostic because the line is saturated) and gas metallic- 
ity (Prochaska et al. 2007) that shows a slope matching 
the local velocity /metallicity relation in dwarf galaxies 
(Dekel & Woo 2003). These correlations may be a con- 
sequence of a relation between the mass of the galaxies 
and their metallicities (Wolfe & Prochaska 1998; Ledoux 
et al. 2006), or an indication that feedback from SF has a 
direct effect on the observed velocity dispersion (Nulsen 
et al. 1998). 

To date all DLAs exhibit metal-line absorption 
(Prochaska et al. 2003), and one expects that they will 
all show significant Mgll equivalent widths. The opposite 
is not true; the majority of Mgll-selected absorbers are 
not DLAs. Unlike traditional low ions, Mgll traces both 
cold and warm neutral material, as well as warm par- 
tially photoionized gas. Recently, Murphy et al. (2007) 
found a correlation between metallicity and the rest- 
frame Mgll equivalent width, suggesting a link between 
kinematics and the metal-enrichment history of the ab- 
sorber. They stressed that the absorption-line kinemat- 
ics should be viewed separately from the host galaxy 
kinematics, in other words, the observed velocity widths 
need not scale directly with the mass of the DLA-hosting 
halos. Bouche et al. (2006) used the cross-correlation be- 
tween 1806 Mgll absorbers and ~ 250, 000 luminous red 
galaxies from the Sloan Digital Sky Survey to find that 
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the absorber halo mass is anti-correlated with the Mg 
II equivalent width, suggesting that at least part of the 
observed velocity dispersion in Mgll regions is produced 
by supernova-driven winds and/or other feedback mecha- 
nisms which could drive cold and warm gas efficiently out 
of lower mass (Mhaio < 10^^-^ M©) galaxies. On the other 
hand, this anti-correlation could also be interpreted as a 
transition from cold to shock-heated gas in more massive 

10^^-^ M0 halos accompanied by a drop in the cumu- 
lative cross-section of neutral clouds found in these halos 
(Tinker & Chen 2007). 

At least for now, observational data are insufficient to 
find a relation between the absorbing halo mass and the 
velocity width vqq of neutral absorption in DLAs. There- 
fore, we should consider at least three distinct physi- 
cal mechanisms which could contribute especially to the 
high-end tail of the DLA velocity distribution: 

1) If cold clouds accreting onto massive galaxies are 
dense enough to survive collisional heating while falling 
into the hot (~ 10^ K) virialized region, they could re- 
tain a sufficiently high neutral fraction with column den- 
sities above 10^^-^ cm~^ (McDonald & Miralda-Escude 
1999). Note that for this mechanism to be viable, the 
accreting gas should be already sufficiently metal-rich, 
as even the most metal-poor DLAs have been found to 
have metallicities Z ^ — 2.8, i.e. it would imply a model 
in which SF was ubiquitous at some time before z = 3 in 
field galaxies of mass 10^ — 10^^ M©. In this scenario a 
large fraction of metals observed in DLAs was produced 
far from massive ( ^ 10^^ M©) halos where we see these 
metals in absorption. On the other hand, it is also possi- 
ble that the low-metallicity accreting material could mix 
efficiently with the more metal-rich environment of the 
massive virialized halo. 

2) Conversely, DLA kinematics could be dominated 
by outflow velocities of supernova-driven winds, i.e. we 
could see feedback directly (Nulsen et al. 1998). This 
hypothesis is indirectly supported by the high detected 
outflow velocities in Lyman-break galaxies (LBGs), the 
higher effective widths VF1526 of the Sill 1526 A transition 
in GRB-DLAs, and the ant i- correlation between the ab- 
sorbing halo mass and the effective line width in Mgll ab- 
sorbers. It would also naturally explain the metallicity- 
velocity correlation, since metals should be produced in 
the same regions which drive the winds. 

3) The apparent inefficiency of in-situ SF in DLAs, 
coupled with their relatively high inferred cooling rates 
suggests a picture in which very compact star-forming 
regions, perhaps associated with LBGs, illuminate much 
more extended neutral gas structures (Wolfe & Chen 
2006). Mailer et al. (2001) used a toy model of DLA 
absorption to show that a large covering factor of the 
cold gas in protogalactic clumps is consistent with obser- 
vations. Moreover, numerical simulations conffim that 
massive halos at 2; = 3 often consist of multiple compact 
galaxies embedded into larger neutral clouds. Combined 
outflows from stellar winds and SNe in these galaxies 
could drive large chunks of surrounding neutral material 
to larger galactic radii where they could pick up a higher 
velocity dispersion from the local galaxy group. 

Although we now have absorption data on ~1000 
DLAs, very few groups have attempted to model DLA 
kinematics in recent years. Both analytical and numer- 
ical models in the literature tend to rely on a velocity 



dispersion put in by hand, in part because of our lack of 
understanding of the physical processes of SF and feed- 
back on sub-galactic scales, and in part due to numerical 
resolution limitations. 

McDonald & Miralda-Escude (1999) used an analytical 
model combining the Press-Schechter formalism with a 
picture of individual spherically symmetric halos consist- 
ing of multiple absorbing clouds to show that the rate of 
energy dissipation corresponding to the velocity disper- 
sion of these clouds necessary to produce the observed 
line profiles far exceeds the rate at which energy can be 
supplied to these clouds by gravitational collapse and 
mergers of halos. They noted that a large fraction of 
DLAs show multiple components in absorption profiles of 
the low-ionization lines, and that this observation can be 
reproduced by a model in which the missing dissipation 
energy comes from supernova explosions. In fact, the 
energy injection rate of 1.8 x 10^^ ergsyr~^/(10^^ Mq) 
reproduces well the fraction of multi-component DLAs 
and the overall absorption line velocity width distribu- 
tion with the median value of ~ 90kms~^. 

However, their model does not provide the physical 
mechanism for transferring the feedback energy to the 
turbulent motions in gas clouds, apart from specifying 
the source of this energy. In addition, it is a highly 
approximate model which assumes spherical exponential 
halos consisting of discrete neutral clouds moving with a 
given velocity dispersion. Moreover, the same feedback 
energy per unit mass is injected into all halos, indepen- 
dently of their mass, an assumption which probably puts 
too much velocity dispersion into small halos. 

On the numerical front, Nagamine et al. (2007) used 
cosmological smoothed particle hydrodynamics (SPH) 
simulations coupled to a phenomenological galactic wind 
model to compute the rate of incidence of DLAs as a 
function of halo mass, galaxy apparent magnitude, and 
impact parameter, for a variable strength of hydrody- 
namical feedback. In their model, gas particles were 
driven out of dense SF regions by hand, by assigning 
a momentum in random directions. The wind mass loss 
rate was assumed to be twice the SF rate, and the wind 
carried a fixed fraction of the SN feedback energy, corre- 
sponding to the fixed wind velocities of 242 km s~^ (weak 
wind) and 484kms~^ (strong wind). Depending on the 
strength of feedback, they found that it evacuated gas 
from predominantly low-mass galaxies and increased the 
cross-section and hence the rate of incidence of more mas- 
sive galaxies. In their weak feedback model ~ 10^'^ Mq 
halos contributed the most to the cross-section, whereas 
with strong feedback the peak shifted to^ 10^^ — 10^^ M© 
halos, depending on numerical resolution, as feedback be- 
came more and more efficient at removing gas from low- 
mass systems. Although, Nagamine et al. (2007) did not 
analyze the velocity width statistics, one would expect 
to see numerous line proffies with vqq > lOOkms"^ in 
their models, due the increased incidence rate of massive 
galaxies. 

The purpose of this paper is twofold. First, we test a 
hypothesis that the observed kinematics is driven primar- 
ily by energy coming from gravitational infall in the pro- 
cess of hierarchical buildup of galaxies. This is essentially 
an extension of the idea put forward by Haehnelt et al. 
(1998) that the observed DLA line proffies are caused by 
a combination of random halo motions, rotation, and in- 
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TABLE 1 
Simulation parameters. 
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fall, coupled to the paradigm of several distinct modes of 
accretion onto growing proto-galaxies (Dekel & Birnboim 
2006). We use the approach developed in our earlier pa- 
per (Razoumov et al. 2006, hereafter Paper I), postpro- 
cessing high resolution adaptive mesh refinement (AMR) 
simulations of galaxy formation with high-angular reso- 
lution radiative transfer of UV ionizing photons. We im- 
proved our algorithm in several ways including a better 
treatment of hydro dynamical heating during the radia- 
tive transfer stage, as well as taking into account SF and 
feedback. We experimented with a number of SF mod- 
els, including the standard four-criterion model of Cen & 
Ostriker (1992) and the two- mode SF model of Sommer- 
Larsen et al. (2003), both of which can be tuned to give 
the observed SF rates in the absence of or with mild feed- 
back. Including strong feedback in our experience tends 
to always suppress an ongoing SF, either unless relevant 
scales in the clumpy interstellar medium (ISM) are re- 
solved, or interaction between the wind and the ambient 
medium has been turned off. 

It is well known that without any special treatment the 
feedback energy is quickly lost away in cooling. Several 
prescriptions have been suggested to alleviate this prob- 
lem; popular algorithms used in particle-based galaxy 
formation models are (1) turning off cooling of gas par- 
ticles in the feedback regions (Thacker & Couchman 
2000; Sommer-Larsen et al. 2003; Stinson et al. 2006) 
and (2) using a subresolution model for the multiphase 
ISM (Springel & Hernquist 2003; Nagamine et al. 2007) 
usually coupled with kinematic feedback in which indi- 
vidual multiphase gas particles are driven out of the star- 
forming regions in random or specified directions. 

Grid-based AMR simulations in which feedback energy 
is simply converted into the gas thermal energy seem 
to be quite effective in limiting SF, although this effi- 
ciency certainly depends on numerical and temporal res- 
olution. In the grid-based simulations presented in Sec- 
tion 2 we do not suppress cooling or use any kinematic 
outflow model, in other words, we enforce a fairly con- 
servative approximation to the role of SF in DLA kine- 
matics, largely limited to conversion of some fraction of 
gas into stars. 

In the second half of this paper (Sec. 3) we examine 
DLA kinematics with TreeSPH models in which radiative 
cooling is suppressed in regions surrounding active SF 
sites. 

2. GRID-BASED AMR MODELS 



To accumulate the absorption line statistics, in order 
to compare models to observations, we need to analyze a 
large cosmological volume, at the same time attempting 
to resolve clumpy gas distribution in individual galaxies. 
To achieve the necessary dynamical range, we use the 
grid-based AMR cosmological structure formation code 
Enzo, running a series of models with box sizes of 4/i~^, 
8/i-\ 16h-^ and 32h-^ comoving Mpc, listed in Table 1. 
All of our models have the base grid size of 128^, with 
up to seven additional levels of refinement everywhere in 
the volume. Each level of refinement features twice the 
spatial and 8 times the mass resolution of the previous 
level. In order to check resolution effects, for each box 
size we also ran one 64^^ base grid model with up to seven 
levels of AMR. In addition, we ran a high-resolution sim- 
ulation HN of the most massive halo in the 32h~^ Mpc 
model HI. For this model we traced the positions of all 
dark matter (DM) particles found at z = 3 inside the 
spheres of radii 2h~^ Mpc and 4h~^ Mpc centered on the 
target halo back to the initial redshift z = 99, and re- 
generated the initial conditions (ICs) with 2X the grid 
and 8X the mass resolution inside the region which con- 
tributed all DM particles to the 4/i~^ Mpc sphere, and 
with 4X the grid resolution and 64X the mass resolution 
inside the region which contributed the DM particles to 
the 2h~^ Mpc sphere. Since the base grid resolution in 
this run is 128^, the effective resolution for the halo is 
512^, with 7 additional levels of AMR, resulting in the 
same physical (200 pc) and mass resolution as the entire 
8h~^ Mpc volume in model Ml. 

Star formation was modeled with discrete stellar par- 
ticles using the prescription of Cen & Ostriker (1992), 
with the overdensity threshold ^cr = 100, the star for- 
mation efficiency egf = 0.1, and the stellar particle mass 
= 10^ Mq. Whenever the conditions for SF were 
met and there was enough gas in a cell to form at least 
one stellar particle, this gas was converted into stars. 
Although a stellar particle was created instantaneously, 
the actual star formation and feedback was modeled 
on the local dynamical timescale. At each timestep, a 
eth = 10~^ fraction of the rest-mass energy of the stars 
produced during that time interval was injected locally 
as thermal energy. 

2.1. UVB radiative transfer 

Our galaxy formation models include a non- 
equilibrium ionization network for hydrogen and helium, 
in the presence of a uniform Madau et al. (1999) UVB. 
The output of these simulations at z = 3 is then postpro- 
cessed iteratively with the UVB radiative transfer code 
FTTE (Razoumov & Cardall 2005) developed for nested 
grids, until an equilibrium position of hydrogen and he- 
lium ionization fronts is found, with the angular reso- 
lution of 192 bins over the 47r sphere, using the same 
UVB as in the hydro dynamical run. This algorithm is 
similar to the one we used in Paper I, but with an im- 
proved treatment of hydrodynamical heating. In Paper 
I we added a location-specific hydrodynamical heating 
term which would keep the temperature of each cell con- 
stant during the iterative radiative transfer stage by ex- 
actly balancing radiative cooling, provided that the ra- 
diation field in that cell stays the same as in the hy- 
drodynamical calculation. Recently we found that due 
to numerical instabilities this approach had a tendency 
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to overestimate the neutral fraction in the shock-heated 
regions of very massive halos where the radiative cool- 
ing timescale is very long. Now we simply do not update 
temperature in shock heated regions with T > Tcr, where 
Tcr = 3 X 10^ K. Although heating by active galactic nu- 
clei via helium photoionization could take the tempera- 
ture above 3 x 10^ K, in practice we found that our results 
are not sensitive to the exact choice of Tcr, as long as it 
stays in the 2x10^ — 5x10^ K range, since shock-heated 
regions have fairly narrow boundaries. 

2.2. Spectrum generation 

To generate and process artificial spectra, we follow 
the procedure outlined in Paper I, with one important 
modification. The principal line width diagnostic used 
throughout this paper is the velocity interval vqq encom- 
passing 90% of the optical depth in the line. This diag- 
nostic is dominated by clouds with large optical depths 
- its more detailed discussion and the comparison to the 
equivalent width can be found, e.g., in Prochaska et al. 
(2007). To measure vqq, we need an unsaturated line 
which also has a sufficient oscillator strength to stand 
out above the noise. In observations such a line is picked 
by hand from a set of low ionization lines associated with 
the absorber. In numerical models we deal with a very 
large number of sightlines to accumulate the statistics 
and therefore need a robust algorithm to pick up the 
suitable line automatically. One solution is to adjust the 
line strength by hand, so that the optical depth in the 
line center is always Tc ~ 2 — 3. Here we simply compute 
vqo prior to applying the noise to the spectrum. This ap- 
proach nevertheless gives us the right measurement since 
f 90 is not sensitive to the choice of the specific line pro- 
vided that all of its components can be seen above the 
noise. 

2.3. Results from AMR models 

Figure 1 provides a quick glance at our results. The 
top panels show the HI column density frequency dis- 
tribution f{N,X) defined such that f{N,X)dNdX is 
the number of DLAs in the intervals [A^, A^ + dN] and 
[A, X -h dX] , where N is the HI column density, and dX 
is the "absorption distance" interval 

dX = ^(l + zfdz. (1) 

The bottom panels of Figure 1 show the line density 
^dla(A), i.e. the number of DLAs per unit absorp- 
tion distance with the velocity width higher than t'go, 
for 64^ + 7 (on the left) and 128^ + 7 (on the right) 
models, respectively. Figure 1 gives the same distribu- 
tions as Figures 4 and 5 in Paper I, with several differ- 
ences. First, our latest models use consistently higher 
grid resolution, across the wider range of volumes from 
4:h~^ to 32/i~^ Mpc. Next, we calculate the neutral hy- 
drogen fraction in shock-heated regions more accurately. 
Finally, we account for conversion of gas into stars, as 
well as feedback, albeit in a conservative form. Only 
model Ml has the same grid and mass resolution as one of 
the runs (C2) in Paper I, however, the above-mentioned 
changes yield slightly lower HI column densities across 
the entire A^hi range. 



At low column densities (10^^-^ — lO^^-^cm"^) we find 
fewer absorbers than observed, with the maximum dis- 
crepancy factor of ^ 1.6 at the DLA threshold. Part 
of this discrepancy may be explained by a Malmquist 
bias in the observational result (O'Meara et al. 2007). 
On the other hand, our models produce consistently 
more high-column density absorbers, missing the turn- 
off at A^Hi ^ lO^^'^cm"^. This behaviour is not sur- 
prising as (1) we do not include the formation of H2 
molecules on dust grains, which is believed to be the 
major sink for neutral hydrogen atoms at higher den- 
sities, and (2) we do not account for gas clumping on 
scales below 100 pc. Apart from the uncertainties of H2 
formation at non-zero metallicities, putting this physics 
in by hand into our current models is not yet practical, 
as we do not resolve the ~ 1 — lOpc scale of individ- 
ual cold clouds in which these molecules form. Note 
that A'hi = 10^^ cm~^ ^ 8 pc~^ is also a critical 
surface density threshold for the formation of the cold 
phase (e.g., Schaye 2004). Our limited numerical reso- 
lution in effect prevents conversion of gas into the cold 
phase which would otherwise trigger gravitational insta- 
bility, i.e. further collapse of densest cores leading to 
much smaller cross-sections. In addition, once the multi- 
phase medium develops, hydro dynamical feedback from 
OB winds and supernovae in such a clumpy medium is 
likely to lead to compression of HI regions into thinner 
shells again producing lower cross-sections. Gas clump- 
ing on these small scales is expected to lower the rate of 
incidence of clouds forming H2 which probably explains 
why the observed transition between the HI and H2 col- 
umn density distributions at z = (Zwaan & Prochaska 
2006) corresponds to fewer absorbers (log / (A", X) ~ —26 
at A'hi = 10^^ cm~^ ) than what is predicted in our cur- 
rent models (log /(AT, A) 24.5 at A^hi = lO^^cm-^). 

We see the same trend found in Paper I that higher 
grid resolution partially alleviates the kinematics prob- 
lem producing a larger fraction of v^q > 70kms~^ ab- 
sorbers. However, the median value of vqq in the resolved 
models Nl, Ml and LI ranges from 39 to 53kms~^, well 
below the observed value of 90kms~^. In model HI the 
median vgo is 64kms~"'^, but this model is clearly not 
resolved as it misses a large fraction of absorbers with 
^haio ^ 10^^ M0 (Fig. 2). Similarly, going to larger com- 
putational volumes at a fixed grid and mass resolution 
(Fig. 3) accounts for progressively more massive halos 
and in smaller boxes tends to produce a larger fraction 
of high- velocity systems. The low median value of vgo is 
further illustrated in the HI column density vs. velocity 
scatter plot (Fig. 4) in which we compare a quasi-random 
sample of 135 observed DLAs to a sample of 135 systems 
picked randomly from each model. 

Let us now look at how halos of different masses con- 
tribute to the overall kinematics. Figure 5 shows the 
vqo velocity widths for those DLA sightlines which pass 
through the virial radii of all DM halos in each model. 
These DLAs do not encompass our entire simulated sam- 
ple, as we see a substantial fraction of A^hi > 2 x 
10^^ cm~^ systems in filaments or intergalactic neutral 
clouds not associated with any particular DM halo (Ta- 
ble 2). This fraction varies with resolution and box size 
but is fairly large in all models, except for the nested ICs 
simulation of the massive halo HN in which there are very 
few intergalactic DLAs. Most of the volume in this model 
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Fig. 1. — HI column density frequency distribution f(N, X) (top panels, see text for details) and line density ^dla(-^) of DLAs with the 
Si II velocity width higher than -usiii vs. -usiii (lower panels) for 4:h~^ Mpc (red lines), 8h~^ Mpc (magenta lines), 16h~^ Mpc (green lines) 
and 32h~^ Mpc (blue lines) models, sd, z = 3. The panels on the left are for the 64^ + 7 (i.e., 7 additional levels of refinement) models, the 
panels on the right are for the higher resolution 128^ +7 models. The dashed lines are fits to the observed distributions at z = 3 (Prochaska 
et al. 2005). 
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Fig. 2. — Connected points show the differential line density 
fiA^/((iX(i log Mhaio) of DLAs as a function of halo mass. Iso- 
lated points on the left (M^alo < 10^ '^Mq) and on the right 
(^halo > 10^^-^ M0) show the cumulative line density dN/dX 
of intergalactic and halo DLAs, respectively. 

is not resolved to have DLA-type absorption, except for 
the 2 Mpc refined region. This region corresponds to a 
higher density peak in the initial density distribution in 
which structures started to form earlier, with less gas left 



TABLE 2 
Fraction of intergalactic DLAs. 

box, Mpc//i low res high res nested 
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over in the diffuse filamentary form, which explains its 
very high fraction of DLAs associated with halos (93%). 

Note that at much higher numerical resolution a lot 
of these intergalactic absorbers would accrete onto low- 
mass halos or form more compact filamentary structures 
with a smaller overall cross-section, so that the values 
given in Table 2 should probably be regarded as upper 
limits to the average fraction of intergalactic DLAs at 
z = 3. 

The scatter in Figure 5 is quite large. Apart from the 
resolution effects, there is a correlation between the ob- 
served median vgo and the circular velocity, although this 
dependence has a fiatter slope than the Vohs ^ O.G^circ 




Fig. 3. — Same as Fig. 1 except that now each panel shows two simulations of same grid and mass resolution but different volumes. 



found by Haehnelt et al. (1998). Moreover, we see that 
this relation might break at higher masses where indi- 
vidual halos consist of multiple components, and unless 
these components are resolved the median of the ob- 
served velocity widths can be as small as ~ 0.2vcirc- 

In the top panel of Figure 6 we plotted the mass- 
weighted gas temperature 

(T)bi„ = J TpdV I J pdV (2) 

averaged over the virial spheres of all halos in a given 
mass bin, as a function of that mass. Only halos con- 
taining at least 100 DM particles are included. It is evi- 
dent that a significant fraction of the gas fails to heat to 
the virial temperature, since its cooling is fairly efficient 
and the infall velocities are relatively small. This result 
was first pointed out by Binney (1977) and more recently 
by Birnboim & Dekel (2003), who showed that in one- 
dimensional models a virial shock does not develop in 
halos below some redshift-dependent critical mass. 

We find that in nested grid models lower-mass halos 
have the mean temperature well above their virial tem- 
perature (Fig. 6) indicating that they are found close 
to the virial radii of more massive halos and are being 
heated by the gas falling into the potential wells of their 
massive neighbors. A fit to all halos with T < Tvir, i.e. 
halos not affected by this ram-pressure heating, produces 
a relationship 

log(T)bin = -3.51 + 0.80 log Mhaio, (3) 

which has a steeper slope than that of the virial temper- 
ature (2/3). This result can be also demonstrated by the 
mass fraction /c of gas at T < 10"^ K (Fig. 6). Due to 
shock heating this fraction decreases with increasing halo 



mass, but in lower-mass (10^ < Mhaio ^ 10^^ M©) sys- 
tems most gas can be actually found in the cold phase. 
As expected, the fraction /c traces the amount of neutral 
gas in halos, as can be seen from the lower panel of Fig- 
ure 6, where we plotted the ratio of the mass- weighted 
neutral fraction 

Mun= f^^^^/ fpdV. (4) 

J nm + nmi / J 

to /c. In order to explain DLA kinematics with large 
velocity widths of low- ionization species, we need a suffi- 
ciently large cold fraction /c of gas in massive halos which 
have a high velocity dispersion. On the other hand, it 
is expected that halos above a certain mass scale may 
not contribute significantly to the observed neutral gas 
absorption due to shock heating of their baryons (Tin- 
ker & Chen 2007). Therefore, our goal is to determine 
precisely the fraction of cold gas inside compact clouds 
and/or accretion filaments which survive shock heating 
while falling into the massive halos. To have a large ve- 
locity dispersion, we need to maximize the number of 
such clouds intersected by any single line of sight going 
through the halo. In our simulations as much as ~ 20% of 
all gas in 10^^ - 10^^-^ Mq halos is at T < 10"^ K (Fig. 6), 
however, this amount of cold material does not automat- 
ically produce the observed velocity dispersion. 

One can argue that the spatial distribution of neu- 
tral clouds inside collapsed halos is equally important 
for probing the velocity space. Could improved numeri- 
cal resolution change this distribution? Without running 
higher resolution models of isolated galaxies and under- 
standing the physics of feedback on the scale of typical 
neutral structures in the ISM of absorber galaxies, it is 
very difficult to extrapolate our results to higher resolu- 
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Fig. 4. — HI column density vs. velocity scatter plot for a sample 
of 135 observed DLAs from the SDSS DLA survey (Prochaska et al. 
2005, top panel) and for 135 randomly picked DLAs from models 
Nl, Ml, LI and HI. 
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Fig. 5. — Median (large symbols) and mean (small symbols) vqq 
velocity widths for DLA sightlines which pass through virial radii 
of DM halos with Mhaio > 10^ M©, vs. Mhaio, in models SI (red 
filled circles). Ml (magenta filled circles), LI (green filled circles), 
HI (blue filled circles) and HN (blue open diamonds). Dotted lines 
show the full range of absorption velocities in each mass bin, sep- 
arately for each model. The thick solid line is the z=3 circular 
velocity. Data points on the left (M^alo < 10^"^ M©) show the ve- 
locity widths for intergalactic DLAs, i.e. absorbers not associated 
with any halo. 
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Fig. 6. — Mass- weighted average gas temperature (top panel), 
fraction of cold (T < 10^ K) gas (center panel), and ratio of the 
mass-weighted average hydrogen neutral fraction to the fraction of 
cold gas (lower panel) inside virial radii of all halos in a given mass 
bin, vs. that mass, for model SI (red filled circles). Ml (magenta 
filled circles), LI (green filled circles), HI (blue filled circles) and 
HN (blue open diamonds). The dotted line in the top panel shows 
the virial temperature at z = 3, while the dashed line shows the 
fit with Eq. 3. Note the data points above this line, indicating gas 
heating in small halos by larger scale dynamics in the vicinity of 
more massive halos. 



tion. Even in our current models resolution effects show 
up in several different ways. First of all, in more mas- 
sive ( > 10^^ M0) halos high mass resolution assists early 
collapse onto a larger number of low-mass DM clumps, 
as opposed to slower accretion via filamentary flows onto 
more massive and more centrally concentrated DM halos 
at somewhat lower redshifts (Fig. 7). Higher grid resolu- 
tion leads to more efficient cooling and in general more 
compact shock heated regions, resulting in a higher sur- 
vival probability of neutral clouds falling on massive viri- 
alized halos. The combined effect of higher mass and grid 
resolution in massive environments is the increased line- 
of-sight velocity dispersion in a larger number of com- 
pact ( < lOkpc) neutral gas clouds, as several of these 
clouds can be normally intersected by a line of sight go- 
ing through such a halo. The resulting spectra feature 
multiple components and the combined velocity widths 
of up to 300kms~^ (Fig. 8), with a median vgo value 
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Fig. 7. — Mass- weighted temperature projection (top panels, in K) and HI column density (lower panels, in cm~^) for the most massive 
halo (8 X 10-*^^ M0) in models HI (refined to the same /max = 7 everywhere in the volume; left panels) and HN (the halo at the center was 
resimulated at 4X the highest grid and 64X the highest mass resolution of HI; right panels). Each panel is 4/i~^ Mpc (comoving) on a side 
and 4h~^ Mpc thick. 



close to 90kms~^ (Fig. 5). 

On the other hand, in lower-mass environments, such 
as isolated < 10^^ halos and filaments, higher nu- 
merical resolution yields a more numerous halo popula- 
tion, and consequently, many more neutral clouds. These 
clouds are sufficiently isolated so that very few of them 
fall on the same line of sight, and the velocity dispersion 
in individual clouds is too small to contribute to the large 
observed velocity tail, independently of resolution. How- 
ever, these lower-mass environments supply the majority 
of low absorbers and are an important contributor to 
the overall column density distribution. 

Ideally, to combine the benefits of models of different 
volumes and resolutions, one would like to convolve the 
incidence rates in Figure 2 with the velocity distributions 
from Figure 5, using the best resolved halos in each mass 
bin in both distributions and adding intergalactic DLAs. 
Note that our 4 — 8 Mpc models represent fairly small 
volumes and are therefore subject to the large cosmic 
variance. Unfortunately, these small volumes are also 
the most expensive to compute, as they start to resolve 
the ISM, and the fraction of the volume which is refined 



adaptively rises steeply. 

3. SPH MODELS WITH EFFICIENT FEEDBACK 

We now turn to high-resolution galaxy formation mod- 
els computed with the TreeSPH code (Sommer-Larsen 
et al. 2003) which, unlike Enzo, uses particle repre- 
sentation of the fluid. SF is similarly modeled with 
a set of discrete star particles which represent a pop- 
ulation of stars born at the same time. Two distinct 
modes of SF depending on the thermal history of the 
star-forming gas are used: (1) early and very efficient 
(egf = 1) SF from gas which has always been cooler than 
Tcr = 3 X lO^K and is found above a critical density 
^H,cr = 0.3 cm~^, and (2) a late mode with much lower 
efficiency (egf = 5 x 10~^) fueled by gas which has been 
heated to temperatures above Tcr and has subsequently 
cooled down to be found in clouds with density above 
^H,cr = 0.01 cm~^. 

For our purposes the most important difference be- 
tween the two codes lies in implementation of feedback. 
In Enzo feedback energy is injected into the thermal en- 
ergy of the ISM which is then allowed to be radiated away 
on a cooling timescale resulting in a very conservative re- 
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Fig. 8. — Distribution of velocity widths along random sightlines 
crossing the most massive halo (8x lO-*^^ Mq) in 32h~^ Mpc models, 
at three different resolutions. 



alization of feedback. In the TreeSPH code feedback is 
modeled in accordance with a given initial mass func- 
tion (IMF), however, cooling is suppressed in the feed- 
back regions for the duration of the starburst, resulting 
in a much more efficient expansion of the heated gas into 
the surrounding medium. This latter approach has been 
shown to preserve a larger fraction of angular momentum 
in the baryonic component as it settles into rotationally 
supported disks, and to lead to more realistic disks at 
z = (Sommer-Larsen et al. 2003). Can these state-of- 
the-art models reproduce the observed DLA kinematics 
at z = 3? 

To address this question, we selected three halos from 
a dark matter-only run of box length 10h~^ Mpc which 
was then resimulated with the TreeSPH code at higher 
resolution in Lagrangian regions enclosing the galaxies, 
with gas (and star) particle masses of 1.1 x 10^ M©. Two 
of these halos (15 and 18) would form massive galaxies 
at z = 0, and the third one (Scl54) would form a clus- 
ter. At z = 2.95 their virial masses are 4.2 x 10^^ M©, 
3.0 X 10^^ Mo and 1.6 x 10^^ M©, respectively. In the 
latter proto-cluster run SF and feedback is modeled with 
a set of discrete star particles assuming a Salpeter IMF. 
One of the proto-galaxies (AY18) was modeled with a 
top-heavy Arimoto-Yoshii IMF, and the other proto- 
galaxy was rerun twice, once with the Salpeter (S15) 
and once with the Arimoto-Yoshii (AY15) IMFs, respec- 
tively. We then took these halos at z = 2.95 and pro- 
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Fig. 9. — HI column density frequency distribution f{N, X) (top 
panels) and line density ^dla(-^) of DLAs with the Si II velocity 
width higher than -usiii vs. -usiii (lower panels) for the two Salpeter 
(S15 and Scl54) and two Arimoto-Yoshii (AY15 and AY18) models 
computed with the TreeSPH code. 

jected them onto a hierarchy of nested grids with inter- 
polated grid resolution of 45 pc for the two proto-galaxies 
and 150 pc for the proto-cluster. We postprocessed these 
datasets with the UVB radiative transfer using the same 
algorithm as for our grid-based simulations (Sec. 2.1), 
and with transfer of stellar ionizing photons using the 
adaptive ray splitting method described in Razoumov 
& Sommer-Larsen (2006). The stellar UV luminosity is 
determined using the population synthesis package Star- 
burst 1999 (Leitherer et al. 1999) with continuous SF 
distributed among all stars younger than 34 Myrs. 

Figure 9 shows the differential HI column density and 
cumulative vgo velocity width distributions for the four 
models. These distributions were computed with the ran- 
dom sightlines passing through the resimulated volumes 
which have the average densities in excess of the mean 
density of the Universe at z = 2.95. Therefore we should 
expect more absorption than from a truly representa- 
tive large cosmological volume, which is evident in the 
upper plot of Fig. 9. However, since these volumes con- 
tain massive environments, they should exhibit a higher 
neutral gas velocity dispersion than a larger cosmologi- 
cal volume. The median vgo velocity width for models 
S15, AY15, AY18, and Scl54 is 32, 34, 31, and 31kms-\ 
respectively. This is a surprising result, as with sup- 
pression of cooling in the feedback regions one would ex- 
pect to see a more efficient conversion of feedback energy 
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into hydro dynamical expansion and, therefore, higher 
outflow velocities. In fact, we see that increasing the 
strength of feedback has some effect on velocities. AY 15 
is the same proto-galaxy as S15, resimulated with the 
top-heavy Arimoto-Yoshii IMF with stronger feedback. 
Figure 9 shows that more energetic feedback removes 
gas from higher column densities and produces a larger 
fraction of high- velocity DLAs, although the median vgo 
changes only slightly from 32 to 34kms~^. 

Note that our kinematics results are independent of 
the assumed UVB. Changing the amplitude and hardness 
of the UVB does not alleviate the lack of high- velocity 
systems, as it affects primarily lower column density ab- 
sorbers, and since there is effectively no correlation be- 
tween the velocity widths and column densities (lower 
panel of Fig. 4), changing the UVB just adds or removes 
both low- and high-velocity DLAs from the simulated 
sample. 

4. DISCUSSION 

We have performed a numerical study of DLA kinemat- 
ics using high-resolution grid-based AMR simulations 
with moderate feedback, as well as particle-based SPH 
models with much more efficient feedback. With grid- 
based models our goal was to test the hypothesis that the 
observed velocity dispersion in DLAs could be explained 
by energy coming from gravitational collapse of cosmic 
structures. The broadest line profiles would then come 
from neutral clouds surviving infall into the massive viri- 
alized ( ^ 10^"^ M©) systems, as a 10"^-^ M© halo at 2: = 3 
has a circular velocity already in excess of 100 km s"""^. 
We used a series of models with a comoving volume 
size ranging from 4/i~^ to 32h~^ Mpc and the maximum 
physical resolution of 100 pc to build a DLA population 
covering halos in the mass range 10^-^ — 10^^ Mq. Al- 
though our models can reproduce well the total incidence 
rate of DLAs at z = 3, our median vqq velocities are of 
order 40 — 50kms~^, well below the observed value of 
90 kms~^. We can point out three possibilities to resolve 
this discrepancy: 

1. More massive environments - The velocity disper- 
sion in individual clouds is a strong function of mass. 
Therefore, one could expect that including more massive 
halos through sampling of the rare density peaks could 
alleviate the velocity problem. On the other hand, all our 
4-16 Mpc high resolution models feature similar velocity 
tails (Fig. 1) even though the mass function in the 16 Mpc 
model LI extends to nearly 50 times more massive halos 
than in the 4 Mpc model Nl. Therefore it seems unlikely 
that the more massive (Mhaio > 8 x 10"^^ M©) and con- 
sequently much rarer halos could make any substantial 
contribution to the overall statistics, even though each 
individual halo could have a velocity dispersion of sev- 
eral hundred kms~^. 

Note that including only massive halos in larger simu- 
lation volumes might produce /(AT, X) at low A^hi which 
is too small (Jedamzik & Prochaska 1998) or too large, 
depending on the technique. In our simulations neglect- 
ing low-mass halos can actually raise f{N^X) at low 
column densities, as all the gas which would otherwise 
be pulled in by low-mass halos stays in extended fila- 
ments which sometimes can produce column densities 
A^Hi > 10^^-^ cm~^. On the other hand, our larger 32 Mpc 
models HO and HI do not have enough grid resolution 



to resolve these filaments, and hence produce small line 
densities ^dla(-^)- 

2. Resolution - Our current simulations are already 
very expensive computationally as we refine adaptively 
everywhere in the volume, and the 128^ base grid models 
have in practice 300^ to 400^ resolution elements. If we 
start from 256^ or a larger base grid, many more halos 
will form drawing in a larger fraction of baryons at earlier 
redshifts, whereas inside individual galaxies we will start 
resolving clumpy interstellar gas. However, without run- 
ning these simulations and sampling a large number of 
galaxies at higher resolution, it is impossible to predict 
reliably the effect on the velocity and column density dis- 
tributions. In addition, higher resolution would allow us 
to better model propagation of supernova- driven winds 
in the clumpy interstellar gas. 

The key question is the relative contribution of lower 
mass ( < 10^^-^ M0) and more massive ( > 10^°-^ M©) 
halos. In our current models these two populations con- 
tribute about equally to the total "galactic" DLA column 
density (Fig. 2). It is possible that at higher resolution 
the relative contribution of the more massive halos might 
grow, as they will feature more numerous spatially sepa- 
rated components falling on the same line of sight, per- 
haps with neutral gas tidal tails and bridges, whereas 
isolated lower mass galaxies could become even more 
compact and pull the remaining cold neutral gas from 
extended filamentary structures decreasing its covering 
factor on the sky. 

Unfortunately, the spatial distribution of absorbing 
clouds cannot be reconstructed from the observed line 
profiles in the velocity space, which would otherwise 
point to the minimum scale needed to resolve muli- 
component DLAs. Lopez et al. (2002) studied an almost 
dust free DLA at z = 2.33 featuring 14 resolved metal 
line components. They found very similar Z/Fe ratios of 
low-ionization species in all 14 components, suggesting 
that these components could originate in clouds physi- 
cally located in a small region of space, perhaps much 
smaller than the resolution limit of our current simula- 
tions. 

3. Local microphysics - The most intriguing possibil- 
ity is that the observed velocity dispersion could arise 
as a result of feedback from SF. More efficient feed- 
back could disrupt the highest column density absorbers 
(A^Hi ^ lO^^-^cm"^) in the lower panel of Fig. 4 cre- 
ating a population of low A^hi systems with high ve- 
locity widths. With such a disruption mechanism put 
in by hand with the galactic wind model, Nagamine 
et al. (2007) saw efficient gas removal from lower- 
mass ( ^ 10^^ M0) halos and increased cross-sections in 
> 10^^ M0 halos. In view of our results, a model in which 
more efficient feedback is obtained without resorting to 
the unphysical assumptions of kinematic winds and/or 
suppressing cooling in feedback regions would seem par- 
ticularly compelling. The key challenge here is to come 
up with a model which describes dynamics of individ- 
ual components of the multiphase ISM separately from 
each other, either resolving them explicitly, or using a 
subresolution model in which different components are 
advected differently on a grid. In these models super- 
nova winds would channel most of their energy into the 
lower density regions between the star- forming clouds. 
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whereas the clouds would be more likely to survive the 
disruptive effect of the winds. 

In Section 3 we used TreeSPH models to test the ef- 
fect that the suppression of cooling has on neutral gas 
kinematics. We found the median velocity widths vqq of 
only 31 to 34 km s~^, even though that these models have 
been shown to expel gas efficiently from the star-forming 
regions at high redshifts and produce realistic galactic 
disks at z = 0. In our view, such low velocity dispersion 
could be explained by a number of factors. 

A. Resolution - When numerical models do not resolve 
density inhomogeneities in the ISM, putting a certain 
amount of supernova feedback into the thermal energy 
of the gas and turning off cooling leads to an efficient 
conversion of this energy into the superwind expansion 
into a nearly uniform medium. However, the mass load- 
ing of the wind in a uniform medium is much larger than 
in a clumpy medium of the same average density, since 
in the latter case the wind will predominantly expand 
into the lower density voids between the clumps creat- 
ing a high velocity outflow which will occasionally sweep 
chunks of neutral material off the edges of denser clouds. 
Clumps of gas entrained in winds are often observed in 
low-redshift outflows Rupke et al. (2005). On the other 
hand, at low resolution winds will displace most gas in 
the galaxy moving it to somewhat higher galactic radii, 
producing "puffy" galaxies with a relatively low velocity 
dispersion. This is exactly what we see in SPH simula- 
tions at 45 pc physical resolution. The expectation is that 
at higher resolution we should see multiphase outflows, 
with higher overall velocities in hot ionized gas, and nu- 
merous dense clouds giving rise to multi-component low- 
ionization absorption lines. Getting such galactic winds 
from first principles is notoriously difficult in numerical 
simulations. In this paper, we argue that the necessary 
condition for obtaining the correct DLA velocity disper- 
sion is a separate dynamical treatment of different com- 
ponents of the multiphase interstellar medium at z = 3. 

In addition to the resolution effects, two other factors 
could have affected our SPH results. Here we just men- 
tion these effects briefly, as clearly their study goes be- 
yond the reach of this paper. 

B. Mean redshift of feedback - A realistic population of 
galaxies at z = can be obtained in cosmological simula- 
tions invoking early {z ^4 — 6) episodes of star formation 
with energetic feedback (Sommer-Larsen et al. 2003). In 
this scenario, the radial distribution of gas and its angu- 
lar momentum arise as a cumulative result of successive 
episodes of star formation driving gas out of the galaxies. 
Other model parameters might result in similar star for- 
mation histories; not all parameters can be constrained 



Binney, J. 1977, ApJ, 215, 483 

Birnboim, Y. &; Dekel, A. 2003, MNRAS, 345, 349 

Bouche, N., Murphy, M. T., Peroux, C, Csabai, L, & Wild, V. 

2006, MNRAS, 371, 495 
Cen, R. &; Ostriker, J. 1992, ApJ, 393, 22 
Dekel, A. &; Birnboim, Y. 2006, MNRAS, 368, 2 
Dekel, A. & Woo, J. 2003, MNRAS, 344, 1131 
Haehnelt, M. G., Steinmetz, M., & Rauch, M. 1998, ApJ, 495, 647 
Jedamzik, K. &; Prochaska, J. X. 1998, MNRAS, 296 
Laursen, P. &; Sommer-Larsen, J. 2007, ApJ, 657, L69 
Ledoux, C, Petitjean, P., Fynbo, J. P. U., M0ller, P., & Srianand, 

R. 2006, A&A, 457, 71 



uniquely from observations. One could speculate that a 
different star formation history, e.g. one featuring ener- 
getic feedback down to redshift 2; ^ 3, would produce 
a higher neutral hydrogen velocity dispersion at such a 
redshift. We note that from the spectra of 2: ~ 3 — 4 
LBGs, outflow velocities of 300 — 400kms~^ are rou- 
tinely inferred (e.g., Pettini et al. 2001; Shapley et al. 
2003). Note, however, also that Laursen & Sommer- 
Larsen (2007) were able to match both the magnitude 
and radial fall-off of Lya surface brightness of typical 
z 3 galaxies, indicating that the spatial distribution of 
neutral hydrogen is correctly predicted by current SPH 
models. 

C. AG N feedback - It is possible that the core of each 
massive galaxy hosts a supermassive black hole which 
might have shown an AGN-type activity at some time 
in the past. Neither grid-based nor particle-based sim- 
ulations presented in this paper include feedback from 
AGNs which could have a profound impact on gas kine- 
matics in host galaxies. 

Feedback from star formation is generally thought to 
play an important role in galaxy evolution. If galactic 
winds give rise to the observed neutral gas kinematics 
in DLAs, the same winds should affect the column den- 
sity distribution. In Paper I and in this paper, we have 
attempted to look at both distributions through large- 
scale, "blind-survey" models, in which we simulate cos- 
mological volumes using aggressive grid refinement for 
all galaxies. However, these models are very expensive 
to compute, especially as we get closer to resolving the 
typical scales of the ISM. Perhaps, a more efficient ap- 
proach is to combine these large-scale numerical surveys 
with simulations of the ISM in isolated galaxies, whether 
or not accounting for mergers and cosmic infall. Fortu- 
nately, there is sufficient amount of observational data 
which can be used to constrain these models. Combin- 
ing traditional QSO-DLA data with a closer look at the 
high-redshift star-forming regions with GRB-DLAs and 
with emission line studies should allow us to learn a lot 
more about young proto-galaxies in coming years. 
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